# Dickstein, Ho, and Mark (2023)
# This script generates fit numbers for counterfactuals.

# * # * # * # * # * # * #
# PRELIMINARIES         #
# * # * # * # * # * # * #
setwd("../../../../library")
source("PreliminariesCode.R")

# Are you using with the fixed markup loop?
markuploop = F
if(markuploop == T){print("You are using the fixed markup loop")}

# Which rating areas?
ravec = 1:7

# Which years?
yrvec = 2016

# What percent of the premium is covered by an employer?
nu = .65

# What is the fixed markup?
if(markuploop == F){
  mrkup <- 0
}

# What is the behavioral share
beshr <- .01

# What alpha-psi cutoff to use:
alphampsicutoff <- .05

# What ACG cutoff to use:
acgHL_boundary <- 1.76

# Where is the counterfactual folder?
if(markuploop == F){
  counterfactual_folder <- paste0(project_folder, "/analysis/counterfactuals")
}

# Which demand specification do you want to use?
if(markuploop == F){
  spectouse <- "lowrisk_simplemh_censor0.2_80sample"
}

# Where is the results file?
if(markuploop == F){
  results_file <- paste0(counterfactual_folder, "/specs/", spectouse, "/output/results.txt")
}

# Set up folder for results
if(markuploop == F){
  dirtosave <- paste0(project_folder, "/analysis/tablesandfigures/release/counterfactuals/", spectouse)
}
dir.create(dirtosave)

# Load parameters from the estimated model or are you supplying your own params?
load_params <- T

# * # * # * # * # * # * # * #
# LOAD DATA AND FUNCTIONS   #
# * # * # * # * # * # * # * #

# Counterfactual Functions
source(paste0(counterfactual_folder, "/library/DA03aRunCounterfactualFunctions.R"))
# Counterfactual Data
source(paste0(counterfactual_folder, "/library/DA03bRunCounterfactualGetdataready.R"))
# Sample Data
left_out_sample <- fread(paste0(project_folder, "/data/explodeddata_leftoutsample.csv"))

# * # * # * # * # * # * # * # * # * #
# RESULTS - COUNTERFACTUAL TYPE ONE #
# * # * # * # * # * # * # * # * # * #

# Base Scenario: No movers
OnlyIndividualMarketCounterfactual <- RunCounterfactual(
  ExplDat = ExplodedData_Ind_small,
  PlanDat = PlanData,
  behaviorShare = beshr,
  maxR = 150,
  subsidytypeInd = "Ratio",
  subsidytypeSG = "Ratio",
  fixed_markup = 0)

# * # * # * # * # * # * # * # * # * # * # * # * # * #
# COMPARING COST AND PLAN CHOICE OF LEFT-OUT SAMPLE #
# * # * # * # * # * # * # * # * # * # * # * # * # * #

# Simulated cost and choice probabilities
foo <- OnlyIndividualMarketCounterfactual$ExplodedDat
is.identified(foo, c("plan_id", "hhid"))
foo[, metal := fcase(
  x_av == 0.6, "bronze", 
  x_av == 0.8, "gold", 
  x_av == 0, "unin", 
  default = "silver"
), ]
dat_est <- foo[, .(
  share.unin.est = sum(share * as.numeric(metal == "unin")), 
  share.bronze.est = sum(share * as.numeric(metal == "bronze")), 
  share.silver.est = sum(share * as.numeric(metal == "silver")), 
  share.gold.est = sum(share * as.numeric(metal == "gold")),
  alpha = mode1(alpha),
  omega = mode1(omega),
  orig_subs_id = mode1(orig_subs_id)
), by = c("hhid")]
dat_est <- unique(dat_est[, -"hhid", with = F])

# Observed cost and choice probabilities
dat_obs <- left_out_sample[year %in% yrvec & choice == 1]
is.identified(dat_obs, "hhid")
dat_obs <- dat_obs[, .(hhid, 
  subscriberid, 
  cost = as.numeric(best_guess_AV > 0) * totpaidmonth_span / 100, # translate to monthly hundreds
  choice_AV = best_guess_AV
), ]  
dat_obs[, cost := ifelse(cost > 89.3707340916662, 89.3707340916662, cost), ]
dat_obs[, metal := fcase(
  choice_AV == 60, "bronze", 
  choice_AV == 80, "gold", 
  choice_AV == 0, "unin", 
  default = "silver"
), ]
prop.table(table(dat_obs$metal))

# Merging together
dat <- merge(
  dat_est, 
  dat_obs, 
  by.x = c("orig_subs_id"),
  by.y = c("subscriberid"),
  suffixes = c(".est", ".obs"),
  all = F)

# Creating estimated cost: 
S <- 100
cost_draw_mat <- matrix(((dat$choice_AV / 100) + ((dat$choice_AV / 100) ^ 2) * dat$omega) * rexp(n = nrow(dat), rate = dat$alpha), nrow(dat), S)
cost_draw_mat <- ifelse(cost_draw_mat > 89.3707340916662, 89.3707340916662, cost_draw_mat)
cost_draw_mat <- ifelse(cost_draw_mat < .2, 0, cost_draw_mat)
dat$cost_est <- apply(cost_draw_mat, 1, mean)

# Creating observered statistics:
col_1 <- with(
  dat,
  c(mean(share.unin.est), 
    mean(share.bronze.est), 
    mean(share.silver.est), 
    mean(share.gold.est), 
    mean(cost_est, na.rm = T), 
    0, 
    mean(ifelse(metal == "bronze", cost_est, NA), na.rm = T), 
    mean(ifelse(metal == "silver", cost_est, NA), na.rm = T), 
    mean(ifelse(metal == "gold", cost_est, NA), na.rm = T)  ))
col_2 <- with(
  dat,
  c(mean(metal == "unin"), 
    mean(metal == "bronze"), 
    mean(metal == "silver"), 
    mean(metal == "gold"), 
    mean(cost, na.rm = T), 
    0, 
    mean(ifelse(metal == "bronze", cost, NA), na.rm = T), 
    mean(ifelse(metal == "silver", cost, NA), na.rm = T), 
    mean(ifelse(metal == "gold", cost, NA), na.rm = T)
  ))
results <- data.frame(
  estimated = col_1, 
  observed = col_2
)
row.names(results) <- c(
  "Share Uninsured", 
  "Share Bronze", 
  "Share Silver", 
  "Share Gold", 
  "Cost - All HHs", 
  "Cost - Uninsured HHs", 
  "Cost - Bronze HHs", 
  "Cost - Silver HHs", 
  "Cost - Gold HHs")
  
# * # * # * #
# SAVING    #
# * # * # * #

fwrite(results, paste0(project_folder, "/analysis/tablesandfigures/release/misc/fit_table.csv"))
